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PREFACE 


This  is  the  first  annual  report  on  Contract  NOQO 14-82-K-0805,  "Nonlinear 
Effects  in  Long  Range  Propagation."  The  study  has  been  organized  under  three  task 
headings  as  follows 

Task  I— Shock  pulse  propagation  in  a  homogeneous  ocean 
Task  II— Nonlinear  propagation  in  a  depth- dependent  ocean 
Task  III— Nonlinear  propagation  in  a  caustic  region 

Considerable  progress  has  been  made  in  all  three  tasks.  In  the  present  report, 
however,  attention  is  focused  on  Task  II.  The  investigation  of  inhomogeneous  ocean 
effects  presented  in  this  report  is  complete,  and  Task  II  is  accordingly  regarded  as 
accomplished. 

A  further  report  is  in  preparation  which  will  summarize  the  substantial 
progress  made  to  date  on  Task  ID.  Also  prior  to  the  final  report  a  paper  is  to  be 
prepared  on  the  Task  I  study,  for  presentation  at  the  10th  International  Symposium 
on  Nonlinear  Acoustics,  Kobe,  Japan  (25-28  July  1984). 


I.  INTRODUCTION 


In  a  real  ocean  whose  properties  are  not  uniform,  but  vary  slowly— on  a 
wavelength  scale— with  position,  the  propagation  of  finite-amplitude  acoustic 
signals  differs  in  three  important  ways  from  that  in  a  uniform  medium. 

(a)  Ray  paths  are  curved  (Fig.  1,  page  5),  and  ray  tube  areas  are  no  longer 

2 

proportional  to  s  .  The  inverse-distance  spreading  law  for  weak  waves  is 
replaced  by 

(S/pc)^  p(t',s)  =  f(t*)  (lossless  medium),  (1) 

where  S  is  the  ray  tube  area,  and 


defines  the  retarded  time  for  outgoing  waves  of  small  amplitude. 

(b)  Cumulative  nonlinear  distortion  of  the  signal  is  described  (as  in  the  case 
of  a  homogeneous  medium)  by  replacing  the  linear  retarded  time  variable 
in  Eq.  (1)  with  a  modified  time  variable  t.  However,  the  nonlinear  part 
of  r  grows  at  a  variable  rate  along  the  ray  path,  which  depends  not  only 
on  S(s)  but  also  on  the  local  fluid  properties  /3(s),  p(s),  c(s). 

(c)  The  rates  of  attenuation  and  dispersion  in  a  real  (lossy)  medium  are 
additionally  dependent  on  position. 

As  a  first  step  towards  a  model  which  includes  all  of  these  effects,  we 
consider  the  following  approximation.  Finite-amplitude  sound  in  a  lossless  ocean  is 
considered  as  propagating  along  the  ray  paths  followed  by  small-signal  waves. 
Self -refraction  of  the  signal  wavefronts  is  thus  neglected.  Furthermore,  the  only 
losses  in  our  model  are  those  which  occur  at  shocks  in  the  waveform. 

Within  this  context,  we  are  able  to  introduce  a  generalized  version  of  the 
reduced  propagation  distance  (x)  to  account  for  effects  (a)  and  (b)  above.  Analytical 


results  have  been  obtained  for  the  reduced  distance  by  assuming  a  horizontally 
stratified  ocean  with  simplified  property  profiles.  To  allow  for  realistic  variations 
of  temperature,  pressure,  and  salinity  with  depth,  a  program  has  been  written  to 
evaluate  x  numerically,  and  results  for  some  typical  profiles  are  presented. 

Finally,  we  note  that  the  importance  of  the  x  generalization  lies  in  the  fact 
that  all  the  predictions  of  weak-shock  theory  for  a  uniform  ocean  become  available 
for  the  nonuniform  case,  merely  through  substitution  of  the  appropriate  x  value. 


n.  LITERATURE  SURVEY 

The  generalization  of  ray  theory  to  allow  for  nonlinear  waveform  steepening  in 

a  slowly  varying  medium  has  been  discussed  in  detail  by  Carlton  and  Blackstock  *  and 
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Pelinovskii,  Petukhov  and  Fridman.  Their  results  are  stated  in  terms  of  a  reduced 
propagation  distance,  denoted  here  by  x.  Physically,  x  represents  the  plane  wave 
propagation  distance,  in  a  medium  of  properties  (po»cq)>  which  yields  the  same 
amount  of  nonlinear  distortion  as  propagation  along  the  actual  ray  path  from  sq  to  s. 
Mathematically,  x  is  given  by  an  integral  expression  along  the  ray  path;  a  simplified 
derivation  is  given  at  the  beginning  of  Section  ID. 

3  4 

Following  an  inconclusive  earlier  paper  by  Fridman,  Petukhov  and  Fridman 
addressed  the  overall  problem  of  describing  blast-wave  parameters  in  a  stratified 
ocean.  They  used  the  weak-shock  model  to  obtain  the  departures,  from  their 
homogeneous  ocean  values,  of  the  peak  pressure  (p  ),  time  constant  ( 0  ),  and  related 
signal  properties.  It  is  important  to  recognize  that  much  of  the  deviation  which 
they  report  from  the  reference  case  is  a  consequence  purely  of  linear  acoustics  (i.e., 
ray  curvature). 

To  explain  this  last  point,  we  note  that  linear  ray  acoustics  predicts  different 

values  of  p  (s) — at  the  same  distance  s--along  different  rays  in  an  inhomogeneous 
* 

ocean.  The  variation  arises  from  the  factor  (S/pc)  in  Eq.  (1),  and  has  nothing  to  do 
with  the  nonlinear  effects  which  are  the  subject  of  the  present  report. 

In  the  present  study  we  avoid  confusion  between  the  linear  and  nonlinear 
effects  of  ocean  property  variations,  by  focusing  specifically  on  the  reduced 
distance.  A  knowledge  of  x  completely  characterizes  the  nonlinear  properties  of  a 
given  propagation  path,  within  the  framework  of  weak-shock  theory.  Consequently 
it  appears  preferable  to  present  information  on  x  for  different  ocean  profiles,  ray 
launch  angles,  etc.,  and  leave  the  user  to  draw  conclusions  appropriate  to  whatever 
particular  initial  waveform  is  of  interest. 
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The  calculation  of  x  in  any  particular  case  requires  a  knowledge  of  the 
nonlinearity  coefficient  (jS),  density  (p),  and  sound  speed  (c)  as  a  function  of 
position.  In  seawater  these  properties  are  related  to  the  local  temperature  (T), 
salinity  (S),  and  pressure  (P).  Convenient  approximations  are  given  in  Refs.  5 
through  9,  and  cover  conditions  encountered  in  most  of  the  earth's  oceans.^  The 
relevant  formulae  needed  for  the  present  study  have  been  assembled  in  Appendix  B. 


.  • 
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III.  ANALYSIS  OF  NONLINEAR  PROPAGATION  IN 
A  SLOWLY  VARYING  OCEAN 

A.  Definition  of  Reduced  Distance 

SOURCE 


TUBE 

FIGURE  1 

DEFINITION  SKETCH  FOR  RAY  TUBE  IN  DEPTH-DEPENDENT  OCEAN 

The  propagation  distance  s  is  measured  from  the  source, 
along  the  curved  ray  path. 

2 

As  indicated  in  the  Introduction,  we  expect  the  quantity  Sp  /pc  (where  p  is  the 
acoustic  pressure)  to  remain  constant  following  a  progressive  wavelet.  The  wavelet 
propagation  speed,  along  the  ray  tube  in  Fig.  1,  is  modified  at  finite  amplitudes  to 

c(s)  +  /3(s)u  “  c  +  ^  P  •  (3) 

Here  /?  is  the  nonlinearity  parameter,  and  u  is  the  particle  velocity. 

A  solution  for  progressive  waves,  to  first  order  accuracy  in  p,  is  therefore  of 
the  form 


■h 


‘f(t'  +  |  p  ds)  , 

pc 


(4) 
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where  the  retarded  time  t'  is  defined  by  Eq.  (2).  Equation  (4)  may  be  reduced  to  the 
standard  form  for  nonlinear  plane  waves  in  a  homogeneous  medium,  by  means  of  the 
substitutions 


and 


(reduced  pressure) 


(reduced  distance)  . 


Here  a  is  a  thermodynamic  property  of  the  medium,  defined  by 


(5) 


(6) 


a  =  Pipe5)'*  . 


(7) 


The  reference  ray  tube  area  SQ  refers  to  a  position  close  enough  to  the  source  that 

the  rays  are  still  straight;  it  will  be  replaced  later  in  terms  of  an  arc  length  (or 

initial  distance)  s  . 

o 


The  solution  in  terms  of  reduced  variables  p,  x  is  thus 


■t 


P  c 
Ko  o 


(8) 


in  which  the  subscript  o  denotes  reference  values,  chosen  for  convenience  to  match 
the  properties  of  the  source  region.  The  reduced  quantities  p,  x  may  be  interpreted 
as  the  plane  wave  pressure  and  propagation  distance  in  the  reference  medium. 


It  is  convenient  for  later  purposes  to  express  the  reduced  distance  x  as 

x=s0«n(iG)  ; 


(9) 
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numerical  results  for  the  dimensionless  factor  G  are  presented  in  the  remainder  of 
this  report.  Note  that  the  function  G  is  equal  to  1  for  a  homogeneous  ocean,  in 
which  case  s  is  simply  the  radial  range  and  Eq.  (9)  reduces  to  the  usual  expression 
for  simple  radial  waves. 

Departures  of  G  from  1  indicate  modification  of  the  effective  propagation 
distance  by  the  combined  effects  of  ray  tube  geometry  and  medium  inhomogeneity. 
In  particular  G  >  1  implies  that  distortion  develops  more  rapidly  than  for  spherical 
waves  in  a  homogeneous  medium,  and  G  <  1  that  distortion  develops  more  slowly. 

B.  Analytical  Formulation  of  x  Integral  in  Terms  of  Ray  Coordinates 

» 

The  ray  tube  area  ratio  in  Eq.  (6)  may  be  expressed  as 

Sin6  cos«0  ,  [i  =  (^)J  (10) 

o 

or  equivalently  as 

^  =  *§  cose  cose„  ,  ]  •  «i) 

O 

Here  0  (see  Fig.  1)  is  the  ray  angle,  and  0q  is  the  ray  launch  angle;  x  is  the 
horizontal  range,  and  xq  (corresponding  to  SQ)  is  the  reference  value  near  the 
source.  These  expressions  assume  a  horizontally  stratified  ocean  and  axisymmetric 
(cylindrical)  spreading.  Equation  (10)  is  the  form  used  by  Petukhov  and  Fridman;** 
Eq.  (11)  is  the  form  used  at  ARL:UT  in  the  MEDUSA  ray  tracing  program.^ 

For  any  given  sound  speed  profile,  the  ray  coordinates  (x,z)  and  the  launch 
angle  derivatives  ( 5  or  t)  may  be  calculated.  Then  if  the  a  profile  is  also  specified, 
Eq.  (6)  for  x  may  be  evaluated  by  integration. 


In  general  the  result  will  not  be  expressible  analytically  in  closed  form. 
However,  some  useful  analytic  results  have  been  obtained  by  assuming  a  power  law 
relation  between  a  and  c,  of  the  form 


(12) 


-n 


These  are  described  below.  For  a  discussion  of  the  validity  of  Eq.  (12),  see 
Appendix  C. 


C.  Explicit  Analytic  Results  for  a  Linear  Sound  Speed  Profile 


By  using  the  known  solution  for  ray  trajectories  in  this  case,  we  find 


-|-  =  (^-)  se c2eQ 
—  s  o 

o  o 


sin  0  -  sin  0 
_ c 

e  -  a 


(13) 


Note  that  s/(0-  0  )  is  a  constant  along  each  ray,  equal  to  the  radius  of  curvature. 
Substituting  in  Eq.  (6),  and  using  Eq.  (12)  along  with 


c  _  COS0 

c  =  cos  0  ’ 

o  o 


(14) 


gives  in  the  limit  of  small  sq 


x  = 


d0 

(sin0-  sin0  ) 
o 


(15) 


Since  the  integral  in  Eq.  (15)  can  be  expressed  in  closed  form  for  n=+l  and  n=0, 

the  function  G  can  be  found  for  these  cases.  In  particular,  the  case  n=l  gives  (using 

12 

formula  2.561(5)  in  Gradshteyn  and  Ryzhik's  Table  of  Integrals  1 


sin0  -  sin0 

G(0o’e)  =  T0“0TcosT 

O 


tan(^- 

tan(-| 


sin0^ 

o 


(16) 


This  result  is  of  interest  since  equation  (12)  with  n=l  gives  a  reasonable  approxi¬ 
mation  to  the  deep  ocean  behavior  of  a  (see  Appendix  C). 


_  • 


The  case  n=0  corresponds  to  the  assumption  of  constant  a  in  the  reduced 
distance  integral,  Eq.  (6).  Although  not  claimed  to  be  realistic,  it  provides  a  useful 
comparison.  The  resulting  expression  for  G  (obtained  using  formula  2.551(3)  in  the 


FIGURE  2 

DISTANCE  MODIFICATION  FACTOR  G  PLOTTED  AGAINST  RAY  ANGLE 
Results  based  on  analytic  solution  for  a  linear  sound  speed  profile. 

Figure  2  shows  G(n=l)  and  G(n=0)  plotted  against  0,  for  an  initial  launch  angle 
of  22.5#(in  the  direction  of  increasing  c).  A  significant  conclusion  may  immediately 
be  drawn:  variation  of  ocean  properties  with  depth  has  an  important  direct  effect 
on  the  reduced  distance  (through  the  a  factor  in  Eq.  (6)),  which  is  comparable  with 
the  effect  of  nonspherical  spreading. 


D.  Asymptotic  Approximations  for  Nonlinear  Sound  Speed  Profiles 


Analytic  results  for  G  may  also  be  obtained  for  more  general  profiles  by 
treating  the  relative  variation  in  sound  speed  as  a  small  quantity.  Thus  we  assume 
the  quantity 


<r  = 


(18) 


to  have  a  numerical  value  much  less  than  1. 


It  is  convenient  to  define  the  sound  speed  profile  in  inverse  form,  by 
expressing  z  as  a  second- degree  polynomial  in  a: 

-  =  a  -  i  btr2  (a,b  constants)  .  (19) 

a  2 

This  allows  a  first-order  (in  a)  departure  from  linearity. 

If  the  x  integral  is  evaluated  using  Eq.  (19),  and  is  consistently  approximated 
to  first-order  accuracy  in  a,  we  obtain  the  following  expression  for  the  factor  G  as 
defined  in  Eq.  (9): 


G  =*  1  -  (n  +  ±  )<r  .  (20) 

Details  of  the  analysis  are  given  in  Appendix  D.  On  the  other  hand,  expanding  the 
available  exact  solutions  for  n=l  (Eq.  16)  and  n=0  (Eq.  (17)  in  powers  of  (0-  0Q)  yields 
the  first-order  approximations 


r  0/2  3 

G(n  =  1)  *  (~)  «  1  -  |  a 

co  2 


r  -1/2  , 

G(n  =  0)«(^-)  *i-i| 


(21) 


(22) 


10 


These  evidently  agree  with  Eq.  (20)  as  special  cases.  Since  neither  the  profile 
curvature  parameter  b,  nor  the  launch  angle  0Q,  appears  in  the  above  results,  a 
useful  conclusion  emerges  from  the  asymptotic  analytical  treatment— namely  that 
the  distance  modification  factor  G  is  principally  a  function  of  the  sound  speed  ratio 
c/cQ,  and  only  secondarily  a  function  of  ray  launch  angle  and  the  precise  shape  of 
the  c(z)  profile.  The  predicted  trend  of  G  with  c/cQ  is  given  by  Eq.  (20)  as 


G 


-(n  +  1/2) 


t 


where  n  is  the  index  in  the  approximate  a  versus  c  relation,  Eq.  (12). 


(23) 


However,  a  condition  attached  to  all  these  asymptotic  results  is  that  the  ray 
angle  must  not  approach  zero  (horizontal  propagation),  since  the  series  expansions  in 
Appendix  D  break  down  at  this  point.  As  a  consequence,  Eq.  (23)  cannot  be  applied 
to  a  ray  which  has  passed  through  a  vertex. 


In  the  remainder  of  this  report,  we  present  exact  analytic  and  also  numerical 
calculations  of  G  for  particular  cases.  By  plotting  G  as  a  function  of  c/cQ  for 
different  launch  angles  and  ocean  profiles,  we  are  able  to  test  the  conclusions 
reached  above. 


E.  Comparison  of  Exact  and  Asymptotic  Predictions  for  the  Linear  Profile  Case 

We  begin  by  testing  Eq.  (23)  against  a  case  for  which  Eq.  (16)  provides  an  exact 
analytic  solution:  namely  the  idealized  linear  sound  speed  profile  with  arc  constant. 
Values  of  the  distance  modification  factor  G,  calculated  from  Eq.  (16),  are  plotted 
in  Fig.  3  against  c/cq.  Both  downward-propagating  rays  (launch  angles  0o=15°, 
22.3*  30*,  43*,  60*,  73*)  and  upward-propagating  rays  (0Q=O*,  -15*,  -30*,  -45*,  -60*, 
-75*)  are  shown. 

The  breakdown  of  the  asymptotic  (c/cQ)"^2  prediction,  given  by  Eq.  (23)  with 
n=l,  is  apparent  for  near-horizontal  rays  as  expected.  Figure  3  shows  that  the 
asymptotic  prediction  is  nevertheless  quite  accurate  for  launch  angles  more  than  30* 
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from  the  horizontal  in  either  direction.  The  curves  for  6o=+60*  and  +75*  are  barely 
distinguishable  from  the  asymptote. 


This  plot,  for  a  linear  sound  speed  profile,  shows  results 

-3/2 

assuming  c  constant.  The  asymptotic  (c/cQ)"  prediction 
is  shown,  together  with  the  exact  analytic  calculation  for 
selected  launch  angles  (as  marked  on  each  curve).  The  rays 
launched  at  15*  and  22.5*  have  been  terminated  at  their 
vertex  points. 
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In  Section  IV  below,  we  describe  a  numerical  scheme  for  evaluating  the  factor 
G  along  any  ray  path,  and  results  are  presented  for  more  realistic  ocean  profiles. 
These  provide  a  more  searching  test  of  the  asymptotic  theory,  and  also  allow  a  study 
of  caustic  formation. 


_  • 


_  • 


.  • 
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IV.  COMPARISON  OF  NUMERICAL  AND  ASYMPTOTIC 
PREDICTIONS  FOR  TYPICAL  OCEAN  PROFILES 


Since  analytic  expressions  for  G,  such  as  Eqs.  (16)  and  (17),  are  available  only 
for  highly  idealized  profiles  of  c(z)  and  a(z),  a  parallel  objective  has  been  the 
development  of  an  efficient  numerical  code  with  which  G  can  be  evaluated 
numerically,  for  any  variation  of  ocean  properties  with  depth. 

The  objective  has  been  met  by  adding  on  to  the  ARLsUT  MEDUSA  program  a 
post -processing  stage,  which  takes  MEDUSA  output  quantities  (such  as  ray  angle, 
path  length,  etc.)  and  performs  a  numerical  integration  along  the  ray  path  to  arrive 
at  G.  The  method  is  described  in  detail  in  Appendix  F.  MEDUSA  is  a  versatile  and 
accurate  ray  tracing  program  developed  at  ARLsUT  by  T.  L.  Foreman  and  described 
in  Ref.  11. 

A.  Accuracy  of  Numerical  Solution 

The  numerical  calculation  scheme  has  been  tested  against  the  known  exact 
solution  for  G,  Eq.  (16),  for  a  linear  c(z)  profile  with  ac=constant.  The  test  case  was 
specified  by 

c  =  cq(1  +  z/a)  ,  a  =  88,800  m;  (24) 

downward-going  rays  were  launched  at  angles  0o=15°,  22.5*,  and  60*  from  depth  z=0, 
and  allowed  to  propagate  to  z  =  2900  m.  Table  I  compares  the  numerical  results  for 
horizontal  range,  x,  and  the  nonlinear  distance  modification  factor,  G,  with  the 
exact  analytical  predictions  (xa,  G&)  for  each  case.  The  quantities  tabulated  are 
x-x_,  in  meters,  and  the  relative  error  |G-G_  1/  I  G_- 1 1  .  Note  that  the  first  of 

a  ■  ■  a  3 

these  indicates  the  accuracy  of  the  basic  ray  tracing  program,  while  the  second  is  a 
measure  of  the  combined  program  accuracy. 

The  first  row  in  Table  I  gives  the  errors  for  the  most  accurate  version  of  the 
numerical  calculation.  Here  the  profiles  c(z)  and  a(z)  are  specified  analytically.  By 
the  end  of  the  propagation  path,  the  accumulated  relative  error  in  G  is  only  of  order 
10"5.  This  is  a  measure  of  the  accuracy  of  the  ray  path  solver,  which  is  further 
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TABLE  I.  ERRORS  IN  COMPUTED  NONLINEAR  RANGE 
FACTOR  G,  AND  HORIZONTAL  RANGE  x 

Propagation  through  a  constant-gradient  layer  at  dif¬ 
ferent  launch  angles  (in  direction  of  increasing  c). 
Layer  thickness  =  2900  m;  sound  speed  ratio  1.033 

(a  =  88  800  m).  The  numbers  N  .  N  denote  the 

c'  a 

number  of  points  used  to  define  the  sound  speed  and  a 
profiles  (see  text).  Results  computed  on  the  ARLsUT 
CYBER  machine  with  32-bit  word  length. 


indicated  by  the  small  errors  in  horizontal  range  (a  few  millimeters,  except  for  the 
shallowest  ray). 

Subsequent  rows  in  Table  I  show  a  progressive  loss  in  accuracy  as  fewer  points 
are  used  to  specify  the  c(z)  and  a(z)  profiles.  Interpolation  of  c  and  ct  values  is 
achieved  by  joining  the  profile  points  with  straight  lines  in  the  (1/c,  z)  plane  and  the 
(a,  z)  plane  respectively.  The  resulting  corners  are  rounded  locally  by  cubic  splines, 
in  order  to  make  the  profiles  continuous  through  the  first  two  derivatives  (i.e.,  no 
jumps  in  profile  slope  or  curvature).  As  a  result  of  the  interpolation  scheme,  even 
though  the  c(z)  input  points  lie  on  a  straight  line,  they  are  not  interpolated  by  a 
straight  line  in  the  program;  and  the  smaller  the  number  (Nc)  of  profile  points,  the 
more  the  interpolated  profile  deviates  from  the  ideal  straight  line  case. 

Table  I  shows  that  reducing  the  number  of  points  in  the  sound  speed  profile 
increases  the  errors  in  x  and  G  (with  the  exception  of  the  60*  ray  path  example, 
where  the  10-point  G  calculation  is  unexpectedly  accurate).  Note  that  the 
increasing  errors  which  appear  as  the  number  Nc  is  reduced  are  not  controlled  by 
the  finite  step  size  in  the  ray  tracing  program.  In  fact  the  step  size  varies 
according  to  the  local  curvature  along  the  ray  path,  being  automatically  selected  to 
meet  an  accuracy  criterion.  As  a  result  the  analytic  c(z)  calculation  requires  fewer 
steps— by  a  factor  of  about  4— than  either  of  the  calculations  for  finite  Nc- 

We  conclude  from  Table  I  that  a  relative  error  in  (G-l)  of  less  than  3  percent 
is  provided  by  the  present  numerical  scheme,  for  ocean  profiles  which  are  specified 
by  at  least  10  points  (and  preferably  by  25  or  more).  The  critical  number  is  N  (the 
number  of  points  defining  the  sound  speed  profile),  since  this  affects  the  accuracy  of 
the  resulting  ray  path. 

The  3  percent  figure  in  fact  represents  a  considerable  achievement,  since 
I  Ga~l  I  is  a  small  number  (varying  in  the  example  from  zero  at  the  beginning  of  the 
propagation  path  to  0.045  at  the  end;  cf.  Eq.  (21)).  It  is  interesting  to  note  that  this 
accuracy  is  degraded  if  care  is  not  taken  in  transferring  MEDUSA  output  data  to  the 
integration  program;  in  order  to  avoid  rounding  errors,  the  data  was  transferred  in 
octal  representation,  with  full  accuracy  retained. 


B.  Specification  of  Ocean  Profiles  for  Numerical  Calculations 


In  order  to  arrive  at  c(z)  and  a(z)  profiles  which  are  consistent  with  real  ocean 
properties,  it  is  convenient  to  begin  with  temperature  and  salinity  profiles  T(z)  and 
S(z).  The  hydrostatic  pressure  P(z)  is  conveniently  given  by  Leroy's  formula,  Ref.  5. 
Appendix  B  summarizes  the  relations  between  these  quantities. 

Since  nonlinear  effects  are  most  likely  to  occur  over  propagation  paths  of 
several  kilometers,  we  are  principally  concerned  with  ocean  properties  in  the  deep 
sound  channel  and  down  to  depths  of  around  10  km.  Conditions  within  the  first 
200  m  below  the  surface,  which  are  much  more  variable,  are  not  important  for  the 
present  study. 

Some  typical  deep  ocean  temperature  and  salinity  profiles,  taken  from 
Refs.  13  and  14,  are  shown  in  Figs.  4  and  5,  respectively.  The  general  pattern  is  a 
rapid  drop  in  temperature  with  increasing  depth  between  about  200  m  and  1000- 
2000  m,  followed  by  a  leveling  off.  Thus  the  temperature  is  almost  constant  at 
greater  depths.  The  salinity  profiles  in  the  North  Atlantic  show  a  trend  similar  to 
the  temperature,  but  there  is  relatively  little  salinity  variation  over  the  whole  water 
column  in  the  North  Pacific  profiles. 

Based  on  these  trends,  "representative"  idealized  profiles  have  been  con¬ 
structed  for  both  T(z)  and  S(z),  made  up  of  2  or  3  straight-line  segments.  The 
profile  points  are  specified  in  Table  II;  they  are  also  marked  on  Figs.  4  and  5.  One 
pair  of  profiles  (open  circles)  represents  the  North  Atlantic,  around  a  latitude  of 
20°  N.  The  other  pair  (solid  circles)  represents  the  North  Pacific,  around  36°  N. 

C.  Sound  Speed  Profiles  and  Computed  Ray  Paths 

The  sound  speed  at  35  specified  depths,  distributed  over  the  water  column 
between  the  surface  and  a  depth  of  10  km,  was  computed  from  the  temperature  and 
salinity  data  in  Table  D.  The  results  are  listed  in  Table  III,  and  the  profiles  are 
plotted  in  Fig.  6.  With  the  tabulated  c(z)  data  as  input,  ray  paths  were  then 
calculated  using  the  MEDUSA  ray  tracing  program.  This  was  done  for  each  of  the 
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TABLE  II.  IDEALIZED  TEMPERATURE  AND  SALINITY  PROFILES 
USED  FOR  NUMERICAL  CALCULATIONS 


Location 

Depth,  z 
(m) 

Temperature,T 

(*C) 

Salinity,  S 
(°/oo) 

N.  Atlantic 

0 

23 

37.0 

(20°  N) 

1,000 

4.5 

35.0 

10,000 

3.5 

34.9 

N.  Pacific 

0 

15 

34.3 

(36*  N) 

600 

4 

34.0 

1,800 

2 

34.6 

10,000 

1 

34.6 

standardized  ocean  profiles,  using  3  different  source  depths;  from  each  source 
depth,  5  or  6  rays  were  launched  at  different  angles.  The  results  are  summarized  in 
Figs.  7,  8,  and  9. 

Figure  7  shows  the  resulting  ray  paths  for  the  North  Atlantic  profile.  For 
source  depths  (zQ)  of  200  m  and  1200  m,  the  shallower  rays  exhibit  a  vertex:  we 
note  that  the  asymptotic  analysis  of  Section  III  gives  no  guidance  beyond  this  point. 
Furthermore,  the  rays  in  Fig.  7(b)  launched  from  near  the  sound  channel  axis 
(zQ  =  1200  m)  at  angles  of  0°  and  7.5*  are  confined  to  the  sound  channel. 

The  channeled  rays  exhibit  caustic  formation  (at  approximately  8  km  and 
29  km  horizontal  range  respectively),  at  which  point  the  concept  of  a  distance 
modification  factor  based  on  nonlinear  ray  acoustics  breaks  down.  Accordingly,  in 
the  results  which  follow  the  calculation  of  G  may  be  terminated  for  one  of  three 
reasons: 

(a)  The  ray  reaches  the  ocean  surface. 
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TABLE  III.  SOUND  SPEED  PROFILES  DERIVED  FROM  THE  IDEALIZED 
TEMPERATURE  AND  SALINITY  PROFILES  OF  TABLE  II 


(a)  North  Atlantic  (b)  North  Pacific  '  • 


// 

Depth  (m) 

Speed  (m/sec) 

// 

Depth  (m) 

Speed  (m/sec) 

1 

0 

1532.260 

1 

0 

1506.533 

2 

10 

1531.942 

2 

10 

1506.106 

3 

20 

1531.622 

3 

20 

1505.677 

4 

30 

1531.299 

4 

30 

1505.244 

5 

40 

1530.973 

5 

40 

1504.809 

6 

50 

1530.644 

6 

50 

1504.370 

7 

70 

1529.979 

7 

70 

1503.483 

S 

90 

1529.302 

8 

90 

1502.584 

9 

110 

1528.614 

9 

110 

1501.672 

10 

150 

1527.203 

10 

150 

1499.811 

11 

200 

1525.375 

11 

200 

1497.413 

12 

250 

1523.472 

12 

250 

1494.937 

13 

300 

1521.494 

13 

300 

1492.382 

14 

400 

1517.307 

14 

400 

1487.035 

15 

500 

1512.807 

15 

500 

1481.373 

16 

600 

1507.986 

16 

600 

1475.460 

17 

700 

1502.842 

17 

700 

1476.423 

IS 

800 

1497.374 

18 

800 

1477.44 7 

19 

1000 

1485.535 

19 

1000 

1479.499 

20 

1200 

1488.710 

20 

1200 

1481.557 

21 

1600 

1495.229 

21 

1600 

1485.690 

22 

2200 

1505.096 

22 

2200 

1494.368 

23 

2800 

1515.071 

23 

2800 

1504.368 

24 

3400 

1525.155 

24 

3400 

1514.486 

25 

4000 

1535.348 

25 

4000 

1524.721 

26 

4600 

1545.645 

26 

4600 

1535.073 

27 

5200 

1556.045 

27 

5200 

1545.539 

28 

5800 

1566.542 

28 

5800 

1556.114 

29 

6400 

1577.130 

29 

6400 

1566.794 

30 

7000 

1587.802 

30 

7000 

1577.572 

31 

7600 

1598.550 

31 

7600 

1588.440 

32 

8200 

1609.362 

32 

8200 

1599.390 

33 

8800 

1620.228 

33 

8800 

1610.410 

34 

9400 

1631.133 

34 

9400 

1621.488 

35 

10000 

1642.065 

35 

10000 

1632.611 
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(Although  ray  tube  areas  could  be  computed  for  the  surface  reflected 
ray,  the  phase  change  on  reflection  makes  it  necessary  to  restart  the 
nonlinear  distortion  calculation  from  this  point.) 

(b)  The  ray  reaches  the  ocean  bottom. 

(Again,  continuation  of  the  G  calculation  beyond  this  point  has  no 
meaning,  except  in  the  special  case  where  the  bottom  reflection 
coefficient  is  equal  to  1.) 

(c)  The  ray  tube  area  approaches  zero. 

(At  the  caustic,  the  ray  tube  area  vanishes  but  G  remains  finite.  An 
asymptotic  analysis  of  the  behavior  of  G  in  this  region  appears  in 
section  IV- F  below.) 

The  program  which  calculates  G  stops  as  soon  as  one  of  these  criteria  is  met. 

Corresponding  ray  paths  for  the  representative  North  Pacific  profile  are  shown 
in  Fig.  8.  The  general  pattern  is  similar  to  that  of  Fig.  7.  Again,  there  are  two 
channeled  rays  shown  in  Fig.  8(b);  the  cycle  lengths  and  ranges  to  caustic  formation 
are  somewhat  larger  than  for  the  corresponding  rays  in  Fig.  7(b). 

In  the  next  three  subsections  we  present  a  discussion  of  the  numerical  G  values 
computed  along  these  ray  paths. 

D.  Numerical  Results  for  the  Distance  Modification  Factor  G:  Variation  with 
Horizontal  Range 

Figures  9  and  10  show  how  the  distance  modification  factor  G  varies  along 
each  ray  path,  for  the  North  Atlantic  and  North  Pacific  cases  respectively.  To 
interpret  these  results,  we  recall  the  definition  of  G.  Given  a  uniform  ocean  which 
matches  the  real  ocean  near  the  source,  and  given  the  same  source  level,  sG  is  the 
propagation  distance  required  to  reproduce  the  signal  distortion  found  at  distance  s 
along  a  ray  path  in  the  real  ocean.  Thus  G  values  greater  than  1  indicate  enhanced 
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nonlinear  distortion,  and  vice  versa.  (Note  that  throughout  this  report,  dissipative 
phenomena— except  at  shocks— are  neglected:  in  reality  a  signal  propagating  along 
a  ray  path  of  increased  length  would  suffer  additional  attenuation  and  dispersion, 
which  are  being  left  out  of  account  for  the  present.*) 

We  first  examine  Fig.  9,  which  presents  results  for  the  North  Atlantic  profile. 
(The  results  for  the  North  Pacific  are  qualitatively  similar- -see  Fig.  10.)  Three 
source  depths  are  used:  (a)  200  m,  (b)  1200  m,  and  (c)  10,000  m,  corresponding  to 
the  ray  paths  in  Figs.  7  and  8.  We  shall  see  that  provided  the  sound  speed  continues 
to  increase  (or  decrease)  monotonically  along  the  ray  path,  the  variation  of  G 
follows  the  trends  predicted  by  the  theory  of  Section  III.  On  the  other  hand,  once 
dc/ds  changes  sign,  G  departs  quite  dramatically  from  the  predictions  of  the 
asymptotic  theory. 

1.  Source  Depth  200  m 

This  case  provides  a  good  illustration  of  the  departures  mentioned  above. 
Figure  9(a)  shows  G  reaching  a  value  of  3.5  as  the  shallowest  ray  ( 0O=7.5°)  reaches 
the  surface.  The  growth  of  G  in  this  case  is  progressive,  starting  from  a  horizontal 
range  of  about  5  km  which  is  the  point  at  which  the  downward- going  ray  crosses  the 
sound  channel  axis.  A  similar  progressive  growth  in  G  may  be  seen  along  the  0Q=15* 
ray,  although  it  is  less  pronounced  (G  reaches  1.8  at  the  surface).  Rays  launched 
more  steeply  downwards  are  cut  short  by  the  bottom,  but  would  presumably  continue 
this  trend  in  an  ocean  of  infinite  depth  (the  present  model  ocean  is  10,000  m  deep). 

Two  principal  observations  may  be  drawn  from  Fig.  9(a).  First,  shallow  angle 
rays  launched  towards  the  sound  channel  axis  can  acquire  sufficient  focusing,  by  the 
time  they  cross  the  axis,  that  quite  large  G  values  accumulate  during  subsequent 
propagation  back  to  the  source  depth.  This  is  a  consequence  of  the  reduction  in  ray 
tube  area  S(s)  relative  to  spherical  spreading. 


*  The  validity  of  the  lossless  fluid  assumption,  for  typical  underwater  shock  pulse 
waveforms  with  various  initial  timescales,  is  being  assessed  (as  a  function  of 
propagation  distance)  as  part  of  a  follow-on  to  the  present  investigation. 
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The  second  observation  relates  to  rays  launched  more  steeply  towards  the 
channel  axis  (e.g.,in  Fig.  9(a),  the  22.5°  ray),  which  penetrate  beyond  the  axis  into 
regions  of  relatively  high  sound  speed.  Somewhat  surprisingly,  in  this  case  the 
initial  focusing  over  the  first  few  hundred  meters  (between  the  source  and  the  sound 
channel  axis)  is  sufficient  to  cause  G  to  go  on  increasing  for  several  thousand  meters 
more,  despite  the  defocusing  effect  of  a  monotonically  increasing  sound  speed  as  the 
ray  propagates  downwards. 

2.  Source  Depth  1200  m 


In  this  case  the  source  is  just  below  the  sound  channel  axis.  Rays 

launched  downwards  therefore  experience  a  progressive  increase  in  sound  speed  until 

a  vertex  is  reached.  According  to  the  asymptotic  theory  of  Section  III,  which 
•  1.5 


predicts  G 


(c/c  / 


one,  and  this  is  observed  in  Fig.  9(b)  (0o=7.5“,  15*,  22.5°). 


,  G  is  expected  to  decrease  slowly  from  its  initial  value  of 

Beyond  the  vertex  the 
asymptotic  theory  offers  no  guidance;  but  the  numerical  results  in  Fig.  9(b)  show  a 
flattening  off,  followed  by  either  a  rapid  increase  as  a  caustic  is  approached  (0q=O#, 
7.5°),  or  a  gentle  upturn  as  the  ray  reaches  the  surface  (©q=  15®,  22.5°). 


The  variation  of  G  near  a  caustic  will  be  discussed  further  in  subsection  F 
below.  Otherwise,  the  main  conclusion  drawn  from  Fig.  9(b) — and  the  similar  results 
for  the  North  Pacific  profile  in  Fig.  10(b)— is  that  G  values  for  sources  close  to  the 
sound  channel  axis  remain  in  the  range  0.85  to  1.0,  caustic  regions  excepted. 

3.  Source  Depth  10,000  m 


The  rays  launched  from  the  deepest  source  location  experience  a 

monotonic  decrease  in  sound  speed  (dc/ds  <  0)  over  the  first  90  percent  or  more  of 

their  path  to  the  surface.  The  asymptotic  theory  of  Section  III  applies  to  this 

region,  and  the  calculated  values  of  G  in  Fig.  9(c)  show  the  predicted  behavior  all 

the  way  up  to  the  sound  channel  axis  (depth  1000  m).  Values  of  G  at  this  depth  vary 

from  1.12  (8=0)  to  1.25  (0  <45°). 
o  o 


if 


Beyond  the  sound  channel  axis,  dc/ds  changes  sign  and  a  one-to-one  relation 
between  G  and  c/cQ  is  no  longer  expected.  Figure  9(c)  shows  that  G  actually 
continues  to  increase  between  the  sound  channel  axis  and  the  surface,  maintaining 
the  trend  established  over  the  earlier  part  of  the  propagation  path. 

Corresponding  plots  of  G  versus  x  for  the  same  three  source  depths  (zq=200, 
1200,  10,000  m)  are  shown  in  Figs.  10  (a,b,c)  for  the  North  Pacific  profile.  The  same 
trends  are  observed  as  in  Fig.  9,  and  the  same  conclusions  apply.  The  main 
differences  of  detail  occur  with  the  caustic  forming  rays  (zq  =  1200  m;  0Q  =  0°, 
7.5*),  and  these  are  discussed  in  subsection  F  which  deals  specifically  with  caustics. 

E.  Numerical  Results  for  the  Distance  Modification  Factor  G;  Relation  between 

G  and  c/c  “ 

_ o 

A  major  conclusion  from  the  theoretical  analysis  in  Section  III  is  that  under 
certain  conditions— principally  when  ray  angles  are  not  too  shallow  and  also  dc/ds 
does  not  change  sign— G  is  a  function  of  c/cq  only.  Some  supporting  evidence  from 
the  present  numerical  study  has  been  discussed  in  the  preceding  paragraphs.  We  now 
test  the  theoretical  conclusion  directly  using  Figs.  11  through  13,  which  plot  G 
versus  c/cQ  for  each  of  the  source  depths  in  turn. 

For  the  shallowest  source  location  (zq  =  200  m),  the  rays  travel  only  a  short 
distance  before  dc/ds  changes  sign  at  the  channel  axis  (z  =  1000  m  for  the  North 
Atlantic,  or  600  m  for  the  North  Pacific  model  profile).  The  theory  is  therefore  not 
severely  tested  in  Fig.  11;  the  plots  terminate  on  the  axis,  before  G  has  changed 
significantly.  Nevertheless  the  collapse  of  results  for  different  launch  angles  (apart 
from  the  shallowest  angles,  ®0=7.5#  and  15°)  is  encouraging. 

The  collapse  is  next  tested  for  rays  launched  from  zq  =  1200  m,  in  Fig.  12. 
The  dc/ds  criterion  terminates  the  plots  at  either  the  channel  axis  (when  0Q=O°),  or 
the  first  vertex  (0o=15°,  22.5°).  The  two  steepest  rays  are  terminated  at  10,000  m 
depth,  giving  G=0.89.  There  is  a  definite  tendency  towards  collapse  on  a  c/cq  basis, 
with  the  main  departures  occurring  as  a  vertex  is  approached. 
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FIGURE  11 

PLOT  OF  G  versus  SOUND  SPEED  RATIO  (log-log  SCALES):  zQ=  200  m 
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FIGURE  12 

PLOT  OF  G  versus  SOUND  SPEED  RATIO  (log-log  SCALES):  zQ  =  1200  m 
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FIGURE  13 

PLOT  OF  G  versus  SOUND  SPEED  RATIO  (log-log  SCALES):  zQ  =  10,000  m 
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Finally,  Fig.  13  shows  the  bottom-launched  ray  results  plotted  in  the  same 
format.  The  plots  in  this  case  all  terminate  at  the  sound  channel  axis.  At  this  point 
the  spread  in  G  values,  across  the  launch  angle  range  -15°to  -80°,  is  only  about 
3  percent  (North  Atlantic  profile)  or  4  percent  (North  Pacific  profile). 

From  these  plots  of  G  versus  c/cQ  we  may  deduce  power  laws  of  the  form 

G  =  <c/corm  ,  (25) 

which  give  a  fairly  accurate  description  of  the  way  the  distance  modification  factor 
varies  along  rays  which  are  not  too  near  the  horizontal.  The  values  of  the  index  m 
listed  in  Table  IV  represent  an  approximate  straight  line  fit  to  the  results  for 
0o=+8O°.  They  range  between  1.05  and  2.3,  as  compared  with  the  value  1.5  arrived 
at  from  an  asymptotic  analysis  in  Section  III  for  an  idealized  ocean  profile. 


TABLE  IV.  INDEX  m  OBTAINED  BY  FITTING  THE  POWER  LAW 
G  =  (c/corm  TO  THE  CURVES  IN  FIGURES  11-13 
(STEEP  RAY  LIMIT) 


Figure  No. 

Index  m 

11  (a) 

1.75 

11  (b) 

1.55 

12(a) 

1.1 

12(b) 

1.05 

13  (a) 

2.3 

13(b) 

2.3 

Depth  range 

200  m  down  to 
sound  channel  axis 

1200  m  to  10,000  m 
(downward  propagation) 

10,000  m  to  channel  axis 
(upward  propagation) 


A  significant  conclusion  from  the  results  presented  in  Figs.  11  through  13  is 
that  an  order  of  magnitude  estimate  of  G  is  possible  on  the  basis  of  c/ cQ  alone, 
provided  the  ray  path  is  such  that  dc/ds  does  not  change  sign.  The  equation 
recommended  for  this  purpose,  namely 

G  =(c/c0)"U5  ,  (26) 


_  • 
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is  a  compromise  which  covers  the  range  of  launch  angles  and  source  depths  in 
Figs.  11  through  13;  actual  m  values  along  the  different  rays  vary  from  0.9  to  2.3. 
Accordingly,  Eq.  (26)  is  expected  to  predict  log  G  (and  hence  the  reduced  distance  x) 
within  a  factor  of  2,  under  typical  ocean  conditions  and  without  requiring  detailed 
information  on  the  ocean  profile  and  ray  path. 

F.  The  Approach  to  a  Caustic 

A  caustic  is  formed  when  neighboring  rays  cross;  the  ray  tube  area  then 
vanishes,  and  according  to  geometric  acoustics  the  intensity  would  be  infinite.  In 
the  region  within  a  wavelength  or  so  of  the  caustic,  this  prediction  is  invalid. 
However,  up  to  this  region  a  marked  increase  in  G  is  expected  as  the  rate  of 
nonlinear  distortion  accelerates  along  the  converging  ray  tube. 

Figure  14  shows  two  neighboring  rays  which  intersect  at  range  xc,  on  the 
caustic.  By  simple  geometry,  we  can  relate  the  area  of  the  converging  wavefront  to 
the  launch  angle  derivative  of  the  ray  angle,  (d0/d0Q),  evaluated  at  the  caustic 
(subscript  c). 


x  -  x 


The  resulting  variation  of  G--which  is  dominated  by  the  ray  tube  area 
convergence— follows  from  the  analysis  given  below. 
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Appendix  F  gives  the  logarithm  of  G  as  an  integral  along  the  ray  path, 


I 


I 


! 


-s 

In  G(s)  =  F(s)  =  j  |a(s')  -  -p-j  ds'  . 


(27) 


The  quantity  A(s)  is  related  to  the  ray  tube  area  S(s),  and  is  defined  by 

1  c  -K  « 

A(s)  =  lim  T-(f")  * 

s-»0  so  Jo  o 
o 


(28) 


It  is  clear  from  Fig.  14  that  as  the  caustic  is  approached  (at  x=x  ),  the  ray 
tube  area  varies  as 


S~  (x  -x)  (x  <  x  )  (29) 

c  c 

and  that  the  other  terms  in  the  integral  above  may  be  regarded  as  constant  in  the 
neighborhood  of  the  caustic.  The  change  in  F(s)  between  x  and  xc  follows  as 


AF 


dx’ 

cosec  ’ 


(30) 


where  K  is  the  constant  of  proportionality  in  the  asymptotic  expression  for  A(s) 
deduced  from  Eqs.  (28)  and  (29): 

A(s)  *  K(x  -  x)-*  ,  (x  <  x  )  .  (31) 

C  (<» 

Although  the  integrand  in  Eq.  (30)  is  singular,  the  integral  is  finite  and  given 
by 


The  constant  K  may  be  deduced  from  Eq.  (28)  and  the  ray  tube  area  relation 

_  xt  cos  9 

=e2  cos6  ’ 

O  0 

o 


I 
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(33) 


if  we  make  use  of  the  asymptotic  result 


(£  =  0  by  definition) 
c 


9F  fey*  I  *  (x  ‘  O 

86o  x  x=x„  c 

c 


Stand 

£0 


X-X, 


“=2*c  *  ( al; )  _<x  -  xc> 


o  c 


(34) 


The  value  of  K  thus  obtained  is 


K  =  jr (cos0  cos ec)>4(‘  3F  )  xc*  ’ 


o  c 


(35) 


which  when  combined  with  Eq.  (32)  above  gives 


«c 
2  — 

o 


(36) 


14 

For  purposes  of  rough  estimation  the  factor  ^acfa0^cc^cc)  maY  replaced 

by  1.  Then  we  find,  on  putting  x  =  0.9  x  ,  that  over  the  final  10  percent  of  the 

♦  ^ 

horizontal  range  up  to  the  caustic,  G  increases  by  a  factor  exp[AF(10%)]  defined 
as  follows: 


where 


G(xc)/G(0.9xc)  =  exp[AF(10%)l  , 


(37) 


(38) 


*  It  is  assumed  that  negligible  ray  bending  occurs  over  the  final  10  percent  of  the 
range,  so  that  the  asymptotic  relation  (36)  may  be  applied.  For  the  relevant  ray 
paths,  see  Figs.  9(b)  and  10(b). 
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Jk. 


The  numerical  examples  discussed  in  this  report  include  four  rays  which  form 
caustics.  These  are  the  rays  launched  at  0Q  values  of  0®  and  7.5*,  from  a  depth  of 
1200  m  (both  ocean  profiles).  For  these  four  rays,  Table  V  lists  values  of  the 
following  parameters:  - 

(a)  The  horizontal  range  xc  at  which  the  caustic  is  formed. 

(b)  The  ray  angle  0  at  the  caustic.  - 

4  W 

(c)  The  derivative  (-d0/d0Q)c,  evaluated  from  the  relation 

(3To)c  =  "x“»2e)c  »»> 

which  follows  from  Eq.  (34). 

(d)  The  factor  exp[AF(10%)]  by  which  G(x)  is  estimated  to  increase 
between  0.9  xc  and  xc>  evaluated  from  Eq.  (38)  above. 


TABLE  V.  RAY  PATH  PARAMETERS  EVALUATED  AT  THE 
CAUSTIC 

Results  for  two  different  ocean  profiles  (as  in 
Table  III).  The  source  is  located  near  the  sound 
channel  axis  (source  depth  =  1200  m). 


®o 

(deg) 

xc 

(km) 

6c 

(deg) 

(-d0/d  0  ) 
o  c 

Profile 

G(xc) 

G(0.9xc) 

0 

7*7605 

-0.01 

3.6747 

N.  Atlantic 

1.39 

7.5 

28*687 

-1.07 

7.2158 

N.  Atlantic 

1.27 

0 

15*4905 

0.15 

5.3047 

N.  Pacific 

1.32 

7.5 

36*730 

-1.14 

6.4757 

N.  Pacific 

1.28 
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The  range  derivative  tx  required  for  the  last  two  table  entries  is  provided  at  each 
point  along  the  ray  path  by  the  MEDUSA  program,  and  values  at  the  caustic  have 
been  deduced  by  linear  interpolation. 

It  is  interesting  to  note,  from  Table  V,  that  (-d0/d0Q)c  is  significantly  greater 
than  one  (of  order  5)  for  the  examples  studied.  The  larger  this  ratio  becomes,  the 
smaller  the  increase  in  G  associated  with  convergence  of  the  ray  tube,  for  a  given 
starting  value  of  x/xc  (assumed  close  to  one). 

The  final  column  in  Table  V  shows  that  the  predicted  increase  in  G,  between 
0.9  xc  and  the  caustic  location  xc,  is  of  order  30  percent.  The  numerically 
calculated  increases,  shown  in  the  plots  of  G  versus  x  (see  Figs.  9(b)  and  10(b) 
discussed  earlier),  are  considerably  less,  of  order  5  to  10  percent.  The  reason  is  that 
the  numerical  plots  are  terminated  at  the  last  step  before  the  caustic  is  reached.  In 
order  to  continue  the  plots  up  to  the  caustic,  it  would  be  necessary  to  modify  the 
MEDUSA  ray  tracing  program  to  locate  the  caustic  automatically. 


V.  CONCLUSIONS 


The  principal  conclusions  from  the  present  study  of  nonlinear  propagation  in  an 
inhomogeneous  ocean  are  as  follows. 

(1)  An  analytical  model  for  finite-amplitude  sound  propagation  in  an  in¬ 
homogeneous  ocean  has  been  set  up,  based  on  ray  theory  combined  with 
the  weak- shock  approximation. 

(2)  Within  the  framework  of  the  model,  separation  has  been  achieved 
between  the  nonlinear  aspects  of  the  problem  and  the  purely  linear 
aspects,  by  defining  a  distance  modification  factor  for  nonlinear  propa¬ 
gation. 


(3) 


The  distance  modification  factor  G  defines  the  equivalent  propagation 
distance  in  a  uniform  ocean, 


sequiv 


=  sG  , 


which  reproduces  the  nonlinear  distortion  found  in  the  inhomogeneous 


case  under  the  same  source  conditions. 


(4)  An  asymptotic  theory,  valid  for  ray  paths  along  which  dc/ds  does  not 
change  sign,  predicts  that  G  is  principally  a  function  of  the  local 
sound  speed  ratio  c/cq. 

(3)  A  numerical  scheme  has  been  developed  for  calculating  G.  Input  data 
for  the  calculation  is  provided  by  a  ray  tracing  program  developed  at 
ARL;UT.  The  scheme  has  been  verified  for  an  idealized  ocean  profile 
with  a  constant  sound  speed  gradient,  by  comparing  the  results  with  an 
exact  analytic  solution  for  G  which  has  been  developed  for  this  case. 

(6)  The  distance  modification  factor  has  been  evaluated  numerically,  for 
two  typical  deep  ocean  profiles  and  several  combinations  of  source 
depth  and  ray  launch  angle.  Up  to  the  point  on  each  ray  path  where 


dc/ds  changes  sign,  the  asymptotic  theory  is  found  to  yield  an  approxi¬ 
mate  collapse  of  all  the  results  in  terms  of  c/cq. 

(7)  Some  of  the  rays  in  the  numerical  study  exhibit  caustic  formation.  The 
increase  in  G  which  occurs,  as  a  ray  approaches  a  caustic,  has  been 
predicted  asymptotically. 

(8)  The  range  limitations  on  weak-shock  theory  are  being  explored  in  a 
separate  study,  with  particular  reference  to  transient  pulses  in  a  real 
ocean  with  chemical  relaxation  effects. 


function  of  s,  Eq.  (F-2) 

profile  constant,  Eq.  (24)  and  Appendices  D,  E;  constant  in  empirical 
Eq.  (C-3)  for  0 

function  of  s,  Eq.  (F- 13) 

profile  constant,  Appendix  D;  constant  in  empirical  Eq.  (C-4)  for  0 
constant,  Eq.  (D-26) 
constant-pressure  specific  heat 
sound  speed 

function  of  s,  Eq.  (F-8);  function  of  (P,T,S),  Eq.  (B-5) 

general  function;  derivative  of  F(s),  Appendix  F 

distance  modification  factor  for  nonlinear  propagation  in  an  inhomo¬ 
geneous  ocean,  Eq.  (9) 

constant  in  asymptotic  A(s)  expression,  Eq.  (31) 
power  law  index,  Eq.  (25);  profile  index,  Eq.  (D- 13) 
number  of  points  used  to  define  numerical  profiles  of  c(z),  a{z) 
power  law  index,  Eq.  (12) 
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hydrostatic  pressure  (kgf/cm  ) 

gauge  pressure 

acoustic  pressure 

reduced  pressure 

ray  tube  area;  salinity  (°/oo) 

distance  along  ray  from  source 

reference  distance  close  to  source 

temperature  (*C) 

time 

retarded  time  for  linear  waves,  Eq.  (2) 
function  of  s,  Eq.  (F-l) 


horizontal  range 
reduced  distance,  Cq.  (6) 
depth  below  surface 
source  depth 

/»(pcV* 

thermal  expansion  coefficient 
nonlinearity  parameter 

launch  angle  derivative  of  z(0Q,x)  on  ray,  Eq.  (1 1) 
ray  angle  (downwards  from  horizontal) 
launch  angle 

cosflo,  Appendix  D;  function  of  (T,S),  Eq.  (B-4) 
launch  angle  derivative  of  x(0Q,z)  on  ray,  Eq.  (10) 
fluid  density 

dimensionless  sound  speed  variation,  Appendix  D 
nonlinear  retarded  time  variable 
latitude  angle 

value  at  source 


value  at  caustic 


APPENDIX  B 

FORMULAS  FOR  HYDROSTATIC  PRESSURE,  DENSITY,  SOUND 
SPEED  AND  NONLINEARITY  PARAMETER  IN  THE  OCEAN. 
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In  the  expressions  below,  P  is  the  hydrostatic  pressure  in  kgf/cnr 


(1  kgf/cnrr  =  0.980665  bar).  Gauge  pressure  relative  to  1  atm  is  denoted  by 

P'  =  P  -  1.033  . 

Also  T  is  the  temperature  in  ®C,  and  S  is  the  salinity  in  °/oo. 


(B—  1) 


Hydrostatic  Pressure 

Leroy's  formula^  lor  hydrostatic  pressure  P  as  a  function  of  depth  z  (in  m)  and 
latitude  p  is  as  follows: 


P  =  1.04  +  0.102506  (1  +  0.00528  sin2p)z  +  2.524E-07  z2  . 


(B-2) 


This  formula  is  supposed  to  be  valid  for  all  oceans  except  the  Black  Sea  and  the 
Baltic. 

Density  (P.T.S) 


The  density  p  (in  kg/m^)  is  given  by  Eckart's  formula^ 


999.97 

p  =  0.698  +  \/F  * 


where 


X  =  1779.5  +  11.25  T  -  0.0745  T*  -  (3.80  +  0.01  T)S  , 


F  =  5890  +  38  T  -  0.375  T  +  3  S  +  0.96784  P  . 


(B-3) 


(B-4) 


(B-5) 


Wilson  and  Bradley'  have  given  an  alternative  formula,  based  on  more  recent  data, 
which  will  be  compared  with  Eq.  (B-3)  in  the  next  phase  of  this  study. 

Sound  Speed  c(P,T,S) 


The  sound  speed  c  (in  m/sec)  is  given  by  Wilson's  formula 


c  =  1449.14  +  Vp  +  VT  +  Vs  +  VSTp 


(B-6) 
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in  which  the  V's  are  polynomials  of  up  to  4th  degree  in  the  indicated  variables. 
Explicit  expressions  are  given  in  Ref.  8. 


Nonlinearity  Parameter  (P,T,S) 

The  nonlinearity  parameter  0  is  given  by  Carlton's  formula^ 

0  =  3.3685  +  9.874E-07  S3 
+  7.799E-06  T3  +  1.027E-03  P* 

-  2.276E-10  P*3  +  5.429E-04  ST 

-  8.641  E-07  STP'  -  1.542E-05  ST2 
+  1.620E-08  T3F  . 


APPENDIX  C 


APPROXIMATE  RELATIONS  FOR  a,  0,  AND  DENSITY  p  IN  SEAWATER 


In  this  Appendix  we  examine  the  validity  of  the  power  law  approximation  for 
a  in  terms  of  sound  speed,  Eq.  (12).  We  also  show  that  the  acoustic  properties 
(a,/9,p)  in  seawater  may  be  quite  accurately  expressed  in  terms  of  just  two 
variables,  sound  speed  and  temperature,  rather  them  the  usual  three  (temperature, 
pressure,and  salinity).  Finally,  the  relation  for  the  nonlinearity  parameter  0  used  by 
Petukhov  and  Fridman  is  shown  to  be  inaccurate,  and  an  alternative  relation  of 
similar  form  is  given. 


Power  Law  Relation  between  a  and  Sound  Speed 


In  Section  III.B,  the  relation 


c  x-n 


-^  =  (— ) 
a  'c 
o  o 


(C-l) 


5  Yi 

was  proposed.  Here  a  denotes  the  product  /3(pc  )"  ,  and  clearly  must  depend  in 
general  on  two  other  variables  (e.g.,  pressure  and  salinity)  besides  the  sound  speed. 
Nevertheless  Eq.  (C-l)  provides  a  first  approximation  to  the  actual  variation  of  a 
with  depth  in  the  ocean,  as  we  shall  show  here. 


Values  of  a  computed  from  the  formulas  given  in  Appendix  B  are  plotted 
against  the  sound  speed  in  Figs.  Cl  through  C3.  Isotherms  are  the  gently  sloping 
curves  marked  T=0,  10,  20  (degrees  C)  at  their  left  ends,  respectively.  Isobars  are 
the  steeply  sloping  curves  marked  P=  1,200,400, ...  at  their  upper  ends,  respec¬ 
tively.  The  reference  value  used  for  normalization  corresponds  to  (10°C,  1  bar, 
35°/oo).  Because  of  the  log-log  scale,  any  power  law  of  the  form  (C-l)  will  appear 
as  a  straight  line  on  the  graph. 

Below  the  first  200  m  in  a  typical  deep  ocean  profile,  e.g., the  North  Pacific 
profile  of  Fig.  4,  the  temperature  falls  relatively  rapidly  from  around  10#C  near  the 
surface  to  2-3*C  at  1-2  km  depth.  Below  2  km  the  temperature  remains  almost 
constant.  The  corresponding  variation  of  with  depth  may  thus  be  deduced  from 
Fig.  Cl,  if  we  assume  a  constant  salinity  of  S=35°/oo.  The  depth  profile  can  be 
traced  on  Fig.  Cl  by  noting  that  the  pressure  (in  bars)  is  approximately  one- tenth  of 
the  depth  (in  meters). 


.  • 


.  • 
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Sound  Speed  —  m/sec 


FIGURE  C  l 

a/«ref  versus  c  (LOG-LOG  SCALES),  FOR  SEAWATER 
WITH  SALINITY  S  -  35°/00 
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Sound  Speed  -  m/sec 


FIGURE  C  2 

a/aref  versus  c  (LOG-LOG  SCALES),  FOR  SEAWATER 
WITH  SALINITY  S-37%o 
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Sound  Speed  -  m/sec 


FIGURE  C-3 

a/araf  versus  c  (LOG-LOG  SCALES),  FOR  SEAWATER 
WITH  SALINITY  S  -  34.5  +  0.1 1  T%0 
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The  resulting  trace  follows  quite  closely  the  line  marked  n  =  l  in  Fig.  Cl.  A 
similar  result  is  obtained  by  tracing  the  North  Atlantic  profiles  of  Fig.  4  onto 
Fig.  C3,  which  incorporates  an  appropriate  temperature-salinity  relation. 

We  conclude  that  the  power  law  Eq.  (C-l),  with  n=l,  provides  a  reasonably 
accurate  model  of  the  depth  dependence  of  a  suitable  for  analytical  studies. 
Under  certain  restricted  conditions  a  different  index  may  be  preferred,  however. 
For  example,  at  0#C  and  depths  less  than  4  km,  n=0  is  more  accurate. 

Acoustic  Properties  as  Functions  of  Sound  Speed  and  Temperature 

In  Fig.  C4,  the  density  of  seawater  is  plotted  against  sound  speed.  The 
same  grid  of  temperature  and  pressure  values  has  been  used  as  in  the  previous 
three  figures,  although  the  constant-pressure  lines  are  not  drawn.  Furthermore 
the  grid  has  been  evaluated  for  four  different  salinity  values  (S=34,  35,  36, 
37°/oo),  although  again  for  clarity  these  are  marked  only  once. 

The  conclusion  drawn  from  Fig.  C4  is  that  the  density  of  seawater  may  be 
determined  quite  accurately  over  a  realistic  salinity  range  from  a  knowledge  of 
the  sound  speed  and  temperature,  without  recourse  to  a  third  parameter  (e.g., 
pressure  or  salinity). 

Figure  C5  shows  that  a  similar  collapse  along  lines  T=const.  occurs  for  the 
nonlinearity  parameter  0,  over  the  same  range  of  salinity.  It  follows  that  a 
knowledge  of  (c,T)  in  the  ocean  is  sufficient  to  estimate  the  important  remaining 
acoustic  properties  0(or  a)  and  p. 

Empirical  Relation  between  Nonlinearity  Parameter  and  Characteristic 
Impedance 

An  expression  for  the  nonlinearity  parameter  of  any  inert  fluid,  in  terms  of 
other  thermodynamic  derivatives,  is 


10 


20 


T=  0 
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FIGURE  C4 

DENSITY  p  versus  c,  FOR  SEAWATER  WITH  SALINITY  S  -  34  to  37°/00 
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Here  a^.  is  the  coefficient  of  thermal  expansion,  and  Cp  is  the  constant- pressure 
specific  heat. 

Guided  by  Eq.  (C-2),  Petukhov  and  Fridman  have  proposed  that  the 
variation  of  p  in  seawater  be  estimated  from 

p  =  1  +  pea  (a  =  const.)  ,  (C-3) 

which  amounts  to  treating  the  bracketed  factor  in  Eq.  (C-2)  as  a  constant. 

Figures  C6  and  C7  show  p  plotted  against  pc  (again  using  the  formulas  in 
Appendix  B),  for  two  different  salinity  values.  Evidently  this  presentation  yields 
quite  an  accurate  single-line  collapse,  for  pressures  and  temperatures  in  the  range 
P  =  1-200  bar,  T=0-20°C.  However,  Eq.  (C-3)  as  it  stands  does  not  fit  the  data  at 
all  closely;  it  predict*;  the  lines  marked  (a)  in  Figs.  C6  and  C7 .  (The  constant  a 
has  been  chosen  to  give  the  correct  value  of  P  at  P  =  1  bar,  T=10#C). 

A  much  better  fit  to  the  data,  for  typical  ocean  conditions,  is  given  by 

P  =  -2  +  peb  (b  =  const.)  .  (C-4) 

Equation  (C-4)  is  represented  by  the  lines  marked  (b)  in  Figs.  C6  and  C7.  A  single 
value  of  the  constant, 

b  =  3*63  •  10'6  m^s/kg  ,  (C-5) 


has  been  used  for  each  salinity. 


FIGURE  C6 

P  versus  CHARACTERISTIC  IMPEDANCE  pc,  FOR  SEAWATER 
WITH  SALINITY  S  =  34%o 
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Impedance  —  kg/sec  mz 


FIGURE  C-7 

0  versus  CHARACTERISTIC  IMPEDANCE  pc,  FOR  SEAWATER 
WITH  SALINITY  S  *  37%o 
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APPENDIX  D 


ASYMPTOTIC  RAY  RELATIONS  FOR  SMALL  TURNING  ANGLES 
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Relations  are  derived  below  for  the  horizontal  range  x,  the  ray  path  length  s, 
the  ray  tube  area  ratio  S/SQ,  and  the  nonlinear  distance  factor  G,  in  terms  of  the 
dimensionless  change  in  sound  speed, 


a 


Here  cq  is  the  sound  speed  at  the  source;  a  basic  assumption  in  the  analysis  is  that  a 
remains  a  small  quantity.  We  also  assume  that  the  gradient  of  c(z)  departs  only 
slightly  from  a  constant  value,  as  explained  below. 


Horizontal  Range 


The  integral 


(D-l) 


may  be  converted  into  an  integral  over  c  (or  <r)  by  introducing  a  relation  between  c 
and  z.  A  convenient  profile  description,  yielding  almost-linear  sound  speed  profiles 
if  the  dimensionless  parameter  b  is  small,  is 


z 

a 


(D-2) 


This  gives  the  differential  relation 


dz  =  a(l  -  bo)d<7  , 


(D-3) 


which  is  used  to  substitute  for  dz  in  Eq.  (D- 1). 


Constancy  of  the  horizontal  phase  speed  with  depth  is  expressed  by 
cos0  =  \(1  +  <r)  ,  (\  =  cos  0Q) 

from  which  it  follows  that 


t,no  1  11  "  *2(1 

tan6  =  \ - rnra — 


(D-4) 


(D-5) 
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We  substitute  this  expression  in  Eq.  (D-l),  along  with  Eq.  D-3),  and  expand  the 

2 

integrand  in  powers  of  the  small  quantity  a.  Retaining  only  or  and  c  terms  gives 


12  2 

x  £*  a  cot  0Q  •  (<r  +  2  (esc  flQ  -  b)<r  ]  ; 


(D-6) 


this  approximation  has  a  relative  accuracy  of  order  a. 


However,  the  asymptotic  result  above  is  not  valid  for  ray  angles  too  close  to 
zero.  The  next-order  (o^)  term,  in  the  brackets  on  the  right  of  Eq.  (D-6),  has  a 
coefficient 


|(cot%o  +  5cot20o)-|  b  ; 


(D-7) 


it  therefore  ceases  to  be  small  compared  with  the  a2  term,  whenever  0  is  of  order 


a  or  smaller. 


Ray  Path  Length 


By  applying  the  profile  assumption  (D-2)  to  the  arc  length  integral 

s-  f  Jfc 
J  sin0  ’ 


(D-8) 


and  making  the  same  approximations  as  for  the  horizontal  range,  we  find 


s  “  a  csc0Q*I  o  +  |  (cot20Q  -  b)a2]  . 


(D-9) 


Again,  this  asymptotic  result  breaks  down  for  ray  angles  which  are  too  shallow.  In 
particular,  it  does  not  hold  through  a  vertex,  since  Eq.  (D-9)  implies  only  one  value 
of  s  exists  at  each  depth. 

More  precisely,  Eqs.  (D-6)  and  (D-9)  imply  that  x  and  s  are  uniquely  deter¬ 
mined  by  the  local  sound  speed;  they  therefore  become  invalid  when  dc/ds  changes 
sign  along  the  ray  path. 
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Ray  Tube  Area  Ratio 


Equation  (9)  gives 


(D-10) 


The  launch  angle  derivative  (£)  of  the  horizontal  range  follows  from  differentiating 
Eq.  (D-6).  Also,  to  first-order  accuracy  in  <r, 


sin0  sin0^  (1  -  cot^0^  •  a)  . 
o  o 

Combining  these  results  yields  the  asymptotic  approximation 


(D-ll) 


at  cot20Q  •  <r2  [  1  +  (csc20Q  -  b)<rJ  . 


(D-12) 


Comparison  of  Asymptotic  and  Exact  Results 

Before  proceeding  to  find  an  expression  for  G,  it  is  instructive  to  check  the 
above  asymptotic  results  against  known  exact  solutions.  Analytic  solutions  are 
available  for  certain  profiles  in  the  family 

■7  m  t 

c  =  c  (1  +^z)  ,  (m  =  const.)1  .  (D-13) 

o  ma 

(a)  Case  m  =  1  (linear  profile): 

x  =  a  sec0  (sin0„  -  sin0)  ;  (D-14a) 

oo 

s  =  a  sec  0Q  ( 0Q  -  0)  ;  (D-14b) 


t  The  integrals  required  for  m  =  1/2, 1/3  are  given  in  Gradshteyn  and  Ryzhik's  Table 
of  Integrals,*2  2*272(3)  and  (7). 
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(D-14c) 


2  2 

r"  sec2 0Q  (sin  0O  -  sine)2  . 

o  o  o 


(b)  Case  m  =  Vi  (parabolic  profile): 


x  =  a  sec  0^  •  I  (sin  0  cos  0  -  sin0  cos0)  +  (0  -  0)1 
o  o  o  o 


(D- 15) 


(c)  Case  m  =  1/3  (cubic  profile): 


x  =  a  sec  0  •  1 3(sin0  -  sin0)  -  (sin  0^  -  sin  0)1 , 
o  o  o 


(D-16) 


In  order  to  check  the  previous  asymptotic  expansions— obtained  using 
Eq.  (D-2)  for  the  sound  speed  profile— against  the  exact  results  given  here,  we  need 
to  relate  the  profile  parameters  b  and  m.  Equation  (D-13)  implies  that 


1+^  =(!♦*> 
ma 


(D-17) 


which  may  be  expanded  in  powers  of  <r  to  yield 


f  =  <t  -  i  (1  -  )tt2  +  0(<73) 


(D-18) 


It  follows  that  to  first-order  relative  accuracy  in  a,  the  profiles  in  Eqs.  (D-2)  and 
(D-13)  are  equivalent  if  we  set 


b  =  1  . 

m 


(D-19) 


Expansion  of  the  exact  results  given  for  cases  (a),  (b)  and  (c)  above— a  rather 
lengthy  process- -demonstrates  agreement  to  first-order  accuracy,  as  expected,  with 
the  asymptotic  expressions  (D-6),  (D-9)  and  (D-12). 

Reduced  Distance 


The  reduced  distance  is  defined  by 


f  <f>  * *'  • 


(D-20) 


sQ  o  o 
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To  convert  the  integration  variable  from  s  to  o,  we  note  that 


ds  = 


dz 

sin0 


f 


then  using  Eqs.  (D-3)  and  (D- 11)  gives 

2 

ds  a  csc0  •  (1  -  boXl  +  cot  0^  •  o)do 
o  o 

2 

«*acsc0  *11  + (cot  0^  -b)o]do. 
o  o 

We  model  the  a(z)  dependence  by  the  relation 

ac11  =  const,  j 


thus 


(1  +  o)~n 


=*  1  -  no, 


(D-21) 


t 


(D-22) 


(D-23) 


(D-24) 


to  first  order  in  o. 

Relations  (D-22),  (D-24),  and  (D-12)  are  now  substituted  in  the  reduced 
distance  integral  (D-20)  to  obtain  the  0(o)  asymptotic  result  (for  sQ— 0) 

x  c*  a  <ro  esc  0q  (In  -Jr  -Co)  ,  (oq  —  0)  ,  (D-25) 

o 

where 


C  =  (n  +  1)  -  |  (esc2  eQ  -  b)  .  (D-26) 

A  more  instructive  version  of  Eq.  (D-25)  is  obtained  by  noting  that  in  the  limit 

o  —O,  the  coefficient  ao„  esc  0^  reduces  to  s^,  while  to  first  order  in  o 
o’  oo  o 


t  Eq.  (D-9)  is  obtained  by  integrating  this  result. 
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-Co  «*in(l  -  Co)  . 
Thus  the  reduced  distance  is  given  in  terms  of  o  by 


(D-27) 


x  »  sQln  I  f  -  (1  -  Co) j  .  (D-28) 

This  result  may  be  compared  with  Eq.  (9)  earlier,  which  defines  the  distance 
modification  factor  G.  The  difference  is  that  Eq.  (D-28)  contains  the  factor  o/oQ  in 
the  argument  of  the  logarithm,  as  compared  with  s/sQ  in  Eq.  (9). 

Equation  (D-28),  with  x  entirely  in  terms  of  o,  would  in  fact  be  appropriate  if 
we  had  chosen  to  follow  Ref.  4  and  examine  departures  of  x  from 

s  In  —  rather  than  s  Jn  — 

o  o - os. 

o  o 

(Note  that  the  two  expressions  coincide  only  for  a  linear  sound  speed  profile  and 
small  turning  angles.)  Petukhov  and  Fridman  in  their  paper**  used  the  first 
expression  as  their  reference;  we  prefer  the  second  because  it  leads  to  a  much 
simpler  asymptotic  result,  as  shown  below. 

The  Distance  Modification  Factor  G 


The  factor  G  is  defined  in  the  present  report  by 

x  =  sin  (^  G)  ;  (D-29) 

so 

its  asymptotic  value  in  the  case  discussed  above  may  be  found  in  terms  of  <r,  by 
equating  arguments  in  (D-28)  and  (D-29)  and  using  the  arc  length  expression  (D-9). 
The  result  is  simply 

G  «  1  -  (n  +  j)o  .  (D-30) 

The  following  points  may  be  noted  in  relation  to  Eq.  (D-30). 

(a)  There  is  no  influence  of  launch  angle  (0Q)  on  the  factor  G  as  defined 
above,  up  to  first  order  in  <r. 


(b)  Likewise  there  is  no  influence  of  profile  curvature  (as  expressed  by  the 
parameter  b),  up  to  first  order. 

(c)  These  features  do  not  apply  if  we  follow  Petukhov  and  Fridman  and 
study  departures  of  x  from  sQln (<y/orQ). 


In  an  earlier  report  on  the  topic  of  nonlinear  sound  propagation  in  a  depth 
dependent  ocean,  Carlton  and  Blackstock*  analyzed  the  case  of  vertical  propagation 
using  a  plane  wave  model.  The  purpose  of  this  Appendix  is  to  point  out  that  even 
the  vertical  propagation  of  sound,  from  a  point  source  in  a  depth  dependent  ocean, 
cannot  be  properly  modeled  without  considering  refraction.  Thus  the  one¬ 
dimensional  analysis  of  Carlton  and  Blackstock,  although  correct  for  the  plane  wave 
model  they  assumed,  should  not  be  applied  to  radiation  from  real  sources  in  the 
ocean,  even  at  large  distances. 

Area  Variation  in  a  Vertical  Ray  Tube 


For  purposes  of  illustration,  we  assume  a  profile  in  which  the  sound  speed 
increases  linearly  with  depth: 


c(z)  =  c  Cl  +  z/a)  . 
o 


(E-l) 


The  horizontal  range  (x)  in  this  situation  is  related  to  the  depth  (z)  and  launch  angle 

(°0)  by 


x  =  a|tan8o  +  [sec20o  -  (z/a  +  l)2  . 


(E-2) 


Figure  El  shows  a  near-vertical  ray  (0Q  close  to  90*),  for  which  the  minus  sign  in 
Eq.  (E-2)  is  appropriate  (the  plus  sign  applies  when  the  ray  has  passed  through  a 
vertex  and  is  traveling  upwards).  Since  cos  0Q  approaches  zero,  it  is  appropriate  to 
write  Eq.  (E-2)  in  asymptotic  form  as 


x  =*  z  cot  0  •  (l  +  z/2a)  +  0(cosJ 0  )  . 
o  o 

Eq.  (E-3)  is  valid  for  all  z  values  provided  z/a  is  of  order  one  or  less. 


(E-3) 


It  follows  that  a  vertical  ray  tube,  formed  by  rotating  the  ray  in  Fig.  E-l 

2 

about  the  z  axis,  has  an  area  (S=t rx  )  which  varies  with  z  as  follows: 

S  .  z  z  2 

r  =  u  +  £)  ,  (z_/z  -*o) .  (E-4) 
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_  1 


Source 


DEFINITION  SKETCH  FOR  NEAR-VERTICAL  RAYS 
Rotation  about  the  z  axis  forms  a  ray  tube  with  its  axis  vertical. 


In  practical  ocean  applications,  z  is  small  compared  with  a— in  other  words  the 
percentage  change  in  sound  speed  along  the  ray  path  is  small— and  Eq.  (E-4)  may  be 
approximated  by 


(E-5) 


in  view  of  Eq.  (E-l). 


Interpretation  of  Eq.  (E-5) 

(a)  In  a  uniform  ocean,  the  c/cQ  factor  equals  unity,  and  Eq.  (E-5)  represents 
spherical  spreading. 

(b)  More  generally,  wavefronts  in  a  stratified  ocean  do  not  spread 
spherically,  even  along  a  vertical  propagation  path. 


(c)  In  calculating  the  reduced  distance  integral 


-  ft  <f f. 

•'0  0 


(E-6) 


along  a  vertical  ray,  it  is  inconsistent  to  allow  for  the  vertical 
variation  of  a  while  neglecting  the  non-spherical  spreading  factor 
c/cQ  indicated  by  Eq.  (E-5)  above. 
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It  is  convenient  to  calculate  the  distance  modification  factor  G  directly, 
rather  than  by  first  calculating  x  from  Eq.  (6)  and  then  inferring  G  from  Eq.  (9).  In 
this  Appendix  we  relate  G  to  a  ray  path  integral  F,  and  then  discuss  the  numerical 
method  used  to  evaluate  F. 

Equation  (6)  for  the  reduced  distance  may  be  written  as 


x  =  I  W(s')  ds'  , 


(F-l) 


where  the  arc  length  weighting  factor  W  is  given  by 


A  '  A  ' 


O  '  O  ‘ 


=  sqA(s)  ,  say. 


(F-2) 


Note  that  A(s)  has  dimensions  of  inverse  length.  It  is  independent  of  the  initial 
distance  sQ  along  the  ray,  in  the  limit  sQ—  0.  By  substituting  for  S/SQ  from  Eq.  (II), 
we  obtain 


o  o 


(F-3) 


The  definition  of  G,  Eq.  (9),  is 


x  =  s  in(G  •  p ) 
°  so 


(F-4) 


which  may  be  differentiated  with  respect  to  s  to  give 


dx  _  so  dG  so 
ds  G  ds  +  s 


Equating  the  right-hand  side  to  sQA(s)— from  Eq.  (F-2)— yields 


(F-5) 


1  dG  ■  i  \  1 

s-s  =  A(s)-?  • 


(F-6) 
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It  follows  that 


G  =  exp  F(s), 


(F-7) 


where 

/•* 

F(s)  =  I  f(sOds'  5  f(s)  =  A(s)  -  .  (F-8) 

*t) 

The  advantage  of  this  formulation,  in  computational  terms,  is  that  the 
integrand  in  Eq.  (F-8)  is  well-behaved  as  s'  tends  to  zero,  unlike  W(s')  in  the  reduced 
distance  integral  (F-l)  which  is  singular.  Moreover,  the  starting  distance  sQ  does 
not  appear  in  Eqs.  (F-7)  and  (F-8). 

A  useful  check  on  the  numerical  values  of  F  obtained  for  small  distances  (or 
small  turning  angles)  follows  from  the  asymptotic  expansion 

f(s)  =  A(s)  -  i  sin  eQ  ,  (s  -sQ)  (F-9) 

which  is  based  on  a  locally  linear  profile  approximation,  as  in  Eq.  (24),  with  otcn 
constant.  The  derivation  is  straightforward  using  Eqs.  (D-9)  and  (D- 12). 
Equation  (F-9)  makes  it  clear  that  f(s)  tends  to  a  constant  value  as  s  tends  to  zero. 


The  numerical  integration  scheme  used  to  evaluate  Eq.  (F-8)  is  described  next. 
Since  the  quantities  needed  to  evaluate  f(s)  are  available  at  each  step  of  the 
MEDUSA  ray  tracing  program,**  it  is  convenient  to  use  the  ray-tracing  steps 
(which  are  of  uneven  length)  as  numerical  integration  steps  in  evaluating  Eq.  (F-8). 
We  use  the  second-order  approximation 


s 


i-1 


f(s)ds  =*  ^  h(f.  +  f.  j)  - 


AhVf'i. 


(F- 10) 


where  h  is  the  step  size  (s.  -  Sj_j)  in  arc  length.  For  the  first  integration  step,  we 


(F-ll) 


1 

f(s)ds  s*  Sjf(Sj) 


since  f(s)  is  asymptotically  constant  as  s  tends  to  zero. 


The  values  of  the  derivative  f *  (s)  required  by  Eq.  (F-10)  are  calculated  as 
follows.  Logarithmic  differentiation  of  Eq.  (F-2)  along  the  ray  path  yields 

f'  (s)  =  A(s)B(s)  +  \  ,  (F- 12) 

sz 

where 

B(s)  =  sin0  +  |  zxzxx  cos30  -  ^  +  -y)  cos0  .  (F-13) 

In  the  above  equation,  x  subscripts  denote  derivatives  with  respect  to  horizontal 
range  (x),  following  the  ray  path.  In  particular, 

zx  =  tan0  .  (F-14) 

The  quantities  needed  to  evaluate  A(s)  and  B(s)  are  provided  by  the  ray  tracing 
program  MEDUSA  at  each  step.  Thus  the  ray  tracing  program  is  run  first,  and  the 
appropriate  output  stored  for  transfer  to  the  integration  program,  which  then 
evaluates  F(sp  at  each  step.  Values  of  G(s-)  finally  follow  from 

G(s.)  =  exp  F(s.)  .  (F-15) 
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